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A hydrodynamic description for an s-component mixture of inelastic, smooth hard disks (two 
dimensions) or spheres (three dimensions) is derived based on the revised Enskog theory for the 
single-particle velocity distribution functions. In this first portion of the two-part series, the macro- 
scopic balance equations for mass, momentum, and energy are derived. Constitutive equations are 
calculated from exact expressions for the fluxes by a Chapman-Enskog expansion carried out to first 
order in spatial gradients, thereby resulting in a Navier-Stokes order theory. Within this context 
of small gradients, the theory is applicable to a wide range of restitution coefficients and densities. 
The resulting integral-differential equations for the zeroth- and first-order approximations of the 
distribution functions are given in exact form. An approximate solution to these equations is re- 
quired for practical purposes in order to cast the constitutive quantities as algebraic functions of 
the macroscopic variables; this task is described in the companion paper. 



I. INTRODUCTION 



Flows of polydisperse particles (mixtures) are ubiquitous in nature and industry alike. Examples of the former in- 
clude pyroclastic flows, landslides, pollutant transport, and planetary rings. Examples of the latter include pneumatic 
conveying of grains, ores, and chemicals; fluidized-bed operation for power production and catalytic cracking; mixing 
of pharmaceutical powders (medication and binder) and poultry feedstock (grains and vitamins). A non- uniform 
particle distribution may be a property of the starting material itself, or it may be intentionally utilized in order to 
improve process performance. For an example of the latter, the addition of fines to a relatively monodisperse material 
has been shown to (i) decrease attrition in high-speed conveying lines (ii) increase conversion in high- velocity, 
fluidized-bed reactors and (iii) improve heat transfer efficiency in a circulating fluidized bed (CFB) combustor 
Polydisperse materials are also known to exhibit counter-intuitive behaviors that have no monodisperse counterpart. 
For example, agitation of polydisperse materials via vibration, free-fall, or flow down an incline leads to segregation 
among unlike particles (de-mixing). Enhancing or suppressing this segregation tendency may be critical to process 
performance, depending on whether the desired outcome is a separated or well-mixed state, respectively. 

In the current effort, attention is restricted to rapid flows, in which particle collisions are assumed to be both 
binary and instantaneous in nature. For monodisperse systems, kinetic-theory-based treatments have been successful 
at predicting not only rapid granular flows (in which the role of the interstitial fluid is assumed negligible) , but have 
also been incorporated into models of high- velocity, gas-solid systems. In particular, kinetic-theory-based descriptions 
are now standard components in both commercial and open-source CFD (computational fluid dynamics) software 
packages for multiphase flows such as Fluent^ and MFIX ( http: / /www.mfix.org/p , respectively. Nonetheless, the 
development and application of kinetic-theory-based descriptions for polydisperse systems is in its infancy relative to 
their monodisperse counterparts, as has been highlighted in several recent review articles and perspectives [3jSS[3]- 
The main challenge associated with the derivation of kinetic-theory-based descriptions for mixtures is the increased 
complexity associated with the additional hydrodynamic flelds and associated transport coefficients, and in particular 
with the accurate evaluation of the collisions integrals. Correspondingly, many of the early contributions resorted 
to assumptions which are only strictly true in the limit of perfectly elastic spheres in a uniform steady state: a 
Maxwellian (single-particle) velocity distribution [1, [^, [l^l or an equipartition of energy [Til . [T^ . [Tsl . [13 . However, the 
presence of a non-MaxweUian velocity distribution in granular flows is well-documented jl5[|l6l[l7[|l8l.[T9l[20l.[2ll.[23|. 
and has been shown to have a signiflcant impact on some transport coefficients f23j. Moreover, a non-equipartition 
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of energy between unlike particles is widely established 24 2a, [2y, [23, [2^, [23, [30, [U, [241 , and has been shown to 
significantly contribute to the driving force for segregation [33| and to lead to a reversal of the segregation direction 
[sj, [3^, [3^1 in certain systems. A more recent theory [s^l involves the lifting of both of these assumptions, except in 
the evaluation of collision integrals involving two unlike particles, in which case a Maxwellian velocity distribution is 
assumed for each particle type. Two current theories exist that do not involve either of these assumptions [ssl. [39|. 
though both are based on the Boltzmann equation and thus are limited to dilute flows. Another key difference 
between existing polydispcrsc theories is the base state used in the Chapman-Enskog (CE) expansion. Some theories 
i, i, M [li, i2, 13^ 14, 37, 39] assume an expansion about a perfectly elastic (molecular equilibrium) base state, 
and thus are restricted to nearly-elastic systems. However, in the CE method the base state must not be chosen a 
priori, but rather it is determined as the solution to the kinetic equation to zeroth order in the gradient expansion. 
This solution is found to correspond to the local homogeneous cooling state (HCS) and was used in Ref. [S^ as the 
reference state to determine the Navier-Stokes transport coefficients of a dilute mixture, without any restriction on 
the level of inelasticity. 

The objective of the current effort is twofold. First, a kinetic-theory-based description for the flow of an s-component 
mixture in d dimensions is derived which (i) incorporates non-Maxwellian and non-equipartition effects, (ii) is appli- 
cable to a wide range of restitution coefhcients, and (iii) is applicable to both dilute and (moderately) dense flows. In 
particular, a CE expansion of the revised Enskog theory for inelastic, hard spheres is carried out for both disks {d — 2) 
and spheres {d = 3) up to the Navier-Stokes order. Second, the derivation of the resulting theory is critically compared 
and contrasted to that of existing theories, in an effort to clearly reveal the implications of various treatments on 
both the governing equations and constitutive relations. For this reason, the derivation is presented in a detailed and 
somewhat pedagogical fashion. This work takes the form of two self-contained, companion papers. In this first paper, 
the results of the exact analysis are given. The follow-on paper details the leading order approximations needed for 
the explicit evaluation of all properties derived here: the distribution functions, the "equations of state" (cooling rate 
and pressure), and the transport coefficients. In addition, the methodology used to obtain these results is critically 
compared there to that of previous theories. 

A confusing issue in the granular community is the context of the Navier-Stokes hydrodynamic equations in freely 
cooling granular gases derived in this paper. The expressions for the Navier-Stokes transport coefficients are not 
limited to weak inelasticity and so the calculations provided here apply even for strong dissipation. The Navier- 
Stokes hydrodynamic equations may or may not be limited with respect to inelasticity, depending on the particular 
states analyzed. The CE method assumes that the relative changes of the hydrodynamic fields over distances of the 
order of the mean free path are small. For ordinary (elastic) gases this can be controlled by the initial or boundary 
conditions. However, in the case of granular fluids the situation is more complicated since in some cases (e.g., steady 
states such as the simple shear flow problem 40]) the boundary conditions imply a relationship between dissipation 
and gradients so that both cannot be chosen independently. In these cases, the Navier-Stokes approximation only holds 
for nearly elastic particles [13] . However, the transport coefficients characterizing the Navier-Stokes hydrodynamic 
equations are nonlinear functions of the coefficients of restitution, regardless the applicability of those equations. 

In spite of the above cautions, the Navier-Stokes approximation is relevant to describe a wide class of granular 
flows. One of them corresponds to spatial perturbations of the HCS for an isolated system. Computer simulations 
have confirmed the accuracy of the Navier-Stokes hydrodynamic equations with their associated transport coefficients 
to quantitatively describe cluster formation [4l|. The same kinetic theory results apply to driven systems as well. 
This is so since the reference state is a local HCS whose parameters change throughout the system to match the 

tsical values in each cell. Another examples of good agreement between theory and simulation |42l ] and experiments 
[4^ include the application of the Navier-Stokes hydrodynamics to describe density/temperature profiles in vertical 
vibrated gases, supersonic flow past at wedge in real experiments psj . and nonequipartition and size segregation in 
agitated granular mixtures [27l . [28l . [4^ . In summary, the Navier-Stokes hydrodynamics with the constitutive equations 
obtained in this paper constitute an important and useful description for many different physical situations, although 
more limited than for elastic gases. 



II. OVERVIEW OF DERIVATION 



The theoretical basis for a hydrodynamic description of molecular gases is most completely established at low 
density using the Boltzmann kinetic equation. There, the CE solution and its prediction of transport coefRcients is 
well-established from both computer simulation and experiment [47| . For a moderately dense gas there is no accurate 
and practical generalization of the Boltzmann equation except for the idealized hard sphere fluid. In that case, the 
Enskog kinetic equation describes the dominant positional corrections to the Boltzmann equation due to excluded 
volume effects of other particles on a colliding pair [13] ■ The neglected velocity correlations are important only at 
much higher densities. The derivation of hydrodynamics and evaluation of transport coefficients based on the Enskog 
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kinetic equation leads to an accurate and unique description of moderately dense gases. The generalization to mixtures 



requires a revision of the original Enskog theory for thermodynamic consistency (revised Enskog theory, or RET) 48 1 , 
and its application to hydrodynamics and mixture transport coefficients was accomplished twenty years ago [49 1. 
As noted above, for granular (dissipative) gases, there remains an open problem of predicting transport properties 
at moderate densities, as occur in current experiments and simulations. This problem is addressed here in its full 
generality using the extension of this revised Enskog theory to inelastic collisions without limits on the number of 
components, densities, temperature, or degree of dissipation. This subsumes all previous analyses for both molecular 
and granular gases, which are recovered in the appropriate limits. 

Due to the extreme length of the derivation, an outline of the steps involved is given here for easy reference. 



• Section III. The starting point of the derivation process is the revised Enskog kinetic equations for mixtures 
of inelastic, hard spheres. These equations for the single-particle, position and velocity distribution functions 
of each species, {/i}, take the form of nonlinear, integro-differential equations, where the integral portion arises 
from the collision operator. 

• Section IV. The macroscopic variables of interest (number density {rti}, etc.) are defined exactly in terms 
of moments of {fi} (e.g., (r) = / d-vfi (r, v,i), where v is the velocity of species i). Thus, the macroscopic 
balance equations can be obtained by appropriate manipulation of the Enskog kinetic equations (e.g., multipli- 
cation by dv followed by integration over the velocity to obtain the species mass balance). At this stage, all 
of the constitutive quantities (cooling rate, stress tensor, conduction, and mass flux) appearing in the macro- 
scopic balances are integral functionals of {fi}, which depend explicitly on space and time only through their 
dependence on {fi}. 

• Section V. In order to obtain a hydrodynamic description (one in which the constitutive quantities are 
determined entirely by the macroscopic or hydrodynamic variables), the concept of a normal solution is intro- 
duced. These are special solutions to the Enskog equations for which the {fi} depend on space and time only 
through an implicit functional dependence on the macroscopic fields (or equivalently as explicit functions of 
these local fields and their gradients at the spatial point of interest) . 

• Section VI. An exact analytical solution for {fi} is not a practical objective in the most general case, and thus 
attention is restricted to states with small spatial gradients. In this case the gradients provide a small parameter, 
allowing a small spatial gradients, or small Knudsen number, expansion (i.e., the CE expansion). The analysis 
is carried out here to first (Navier-Stokes) order: fi ~ ff^^ + f'f'\ where ff'' is the zeroth order solution and 
f^^^ is the first-order correction (zero- and first-order in gradients, respectively). The kinetic equations then 
become integral-differential equations for the determination of ff*^ and f\^\ 

• Section VII. Correspondingly, the constitutive equations are identified as functions of the hydrodynamic 
variables and their gradients through their dependence on {fi}, with coefficients of the gradients defining the 
transport coefficients. Hence, all equations of state (pressure and reference state cooling rate) and all transport 
coefficients, which are integrals involving ff^\ inherit this dependence on the hydrodynamic variables and their 
gradients. The coefficients are determined from solutions to the integral equations. 

This completes the derivation reported in this manuscript. Up until this point, the results are exact for Navier- 
Stokes order hydrodynamics (first order in spatial gradients) of the RET. This determines the form of the Navier-Stokes 
hydrodynamics, but more explicit dependence of the transport coefficients on the macroscopic variables requires a 
corresponding explicit solution to the integral equations for f^^^ and /f^"*. One approximate method, known to be 
accurate for ordinary fiuids, is detailed in the follow-on paper [50|, resulting in constitutive quantities that are algebraic 
functions of the macroscopic variables. 



III. REVISED ENSKOG KINETIC THEORY 



The system considered is a mixture of {Ni} smooth hard disks {d — 2) or spheres (d — 3) of masses {rrii} and 
diameters {ui}, where the subscript i labels one of the s mechanically different species and d is the dimension. In 
general, collisions among all pairs are inelastic and are characterized by independent constant normal restitution 
coefficients {aij = ciji}, where a^- is the restitution coefficient for collisions between particles of species i and j, 
< Uij < 1. The macroscopic (or hydrodynamic) properties of interest (number densities, flow velocity, and energy 
density) are determined from the single particle position and velocity distribution functions /i(ri, vi; t), for i = 1, ..s, 
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where /^(ri, Vi; /:)dridvi is proportional to the probabihty to find a particle of species i in the position and velocity 
element dridvi at time t. The fundamental description of any system is based on the probability density for all 
constituent particles and the Liouville equation for its time evolution; this is equivalent to solving the collective 
equations of motion for all particles in the system and becomes computationally prohibitive for a large number 
of particles. However, for the macroscopic fields only the reduced distribution functions {fi}, obtained from the 
integration of the probability density over all except one particle's position and velocity for each of the species, are 
required for calculation of the macroscopic properties. The equations for these reduced distribution functions resulting 
from the partial integrations of the Liouville equation, give rise to the BBGKY hierarchy equations. The first level of 
this hierarchy gives the time dependence of {fi} jsil. [53| 

{dt + VI • Vri+mriF,(ri) • VvJ /^(ri, vi;i) = ^(ri, vi;i), (3.1) 

where 

a(ri,vi;t) = iZK^' J / d^e(^-gi2)(^-gi2) (3.2) 

X ("yVy(ri,v'/,ri - cry,v^';i) - /y(i-i,vi,ri + cry,V2;0) ■ 

The left sides of these equations describe changes in the distribution functions due to motion in the presence of 
external conservative forces Fi(ri). The right side describes changes due to collisions among the particles. The 
function /ij(ri, Vi, r2, V2; t)(iri(ivi(ir2(iv2 is proportional to the joint probability of finding a particle of species i in 
dridvi and one of species j in dr2d-V2- The position r2 in these functions appears only for r2 = ri ± crij, where 
(Tij — aaij and aij = (ui + Gj) /2; this means that the two particles are at contact. The vector is a unit vector 
directed along the line of centers from the sphere of species j to that of species i at contact and the integration 
da is over a solid angle for this contact sphere. The Heaviside step function assures that the relative velocities 
gi2 = Vi — V2 are such that a collision takes place, and the "restituting" (pre-collisional) velocities v'/ and V2 are 
related to the post-collisional velocities by 

v'/ = vi - (1 + a^r.i) (ct • gi2)CT, v^' = V2 + (1 + a^-") [d- ■ §12)^ (3.3) 

where /i^j = mi/ {nii + rrij). It is convenient for the discussion here to write that equation in a more symbolic form 
by introducing the notation 

^(v'i',vD = 6r.ix(vi,V2), (3.4) 

so that b^j^ is a general substitution operator that changes the argument of a function to its precoUision velocities 
given by p.3|) . Then, changing variables a —a in the second term on the right side of (|3.1|) and noting that 
• gi2 = • gi2 gives the equivalent form [sil, (Ell 

(9t+vi.Vr,+mriF,(ri). VvJ/,(ri,vi;t) = / f da{aT^'bT^' + 1)0 ■ ^^2) 

j=i J J 

xe(-CT • gi2)/y(ri, vi,ri - (Tij,V2;t). (3.5) 

This demonstrates that the two particle distributions fij appear only on the contact hemisphere given by 0(— • gi2), 
correponding to particles that are directed toward each other and hence have a change in their velocities. 

Equation (|3.ip becomes a kinetic theory (i.e., closed equations for the set of fi) only after specifying fij on the 
right side as a functional of the set of fi (the alternative of making approximations at higher levels of the BBGKY 
hierarchy has not been productive in general for molecular gases). As indicated above, this is required for fij only 
when the particles are at contact and on that hemisphere for which the relative velocities are directed toward each 
other. In this restricted context, the Enskog kinetic theory results from a neglect of velocity correlations, i.e. the 
Enskog approximation 

/y(ri,vi,r2,V2;t) (ri,r2 | {n,}) /,(ri, Vi; i)/j(r2, V2; i). (3.6) 

Spatial correlations arising from volume exclusion effects are retained through the factor Xij (ri,r2 I {'^i})- In the 
special case of a uniform system, it is simply related to the nonequilibriurn pair correlation function gij (|ri — r2| ; {fii}) 
(probability density to find a particle of species i at ri and j at V2) by [53| 

1 + a ■ 

"^3 
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This relationship is proved in Appendix \^ and provides some partial interpretation for Xij ■ It is important to note 
that these correlation functions are functionals of the actual species densities {n-i} (defined below in Eq. (j4.ip ). This 
functional dependence is what distinguishes the RET from the original "standard" Enskog theory (SET), where the 
gij are functions of the species densities at the single position of interest, ri. Some partial justification for the 
approximation (|3.6[) for ordinary atomic fluids is given in Appendix |^ where it is known to provide accurate results 
for moderately dense gases, and reasonable estimates even for dense gases. Its use for granular gases is justified largely 
from expectations based on these results for ordinary fluids. 

Substitution of the Enskog approximation p.6p into the exact first level hierarchy equations p.ip defines the RET 
for the distribution functions {/;} 

s 

(at+vi.V+m-iF,(ri).VvJ/.(ri,vi;t) = ^Jy[ri,Vi 1/(0] • (3.8) 

The collision operators {Jij [ri,vi | f{t)]} are given by 

Jy[ri,vi I /(t)] = afr^ J dv2 J d^e(^ • gi2)(^ • gia) 

X ['^^j^X^] (ri,ri - (T,j I {n,}) /,(ri, v'/; i)/j(ri - cr,y,v^';t) 

-Xjj (i"i,ri +crjj I {n,})/,(ri,vi;i)/,(ri +(Ty,v2;t)] . (3.9) 

The corresponding Boltzmann equations for a dilute mixture follow from this result since Xij (^i, ri — ""ij | {n-i}) 1 
at low density. Furthermore, on length scales of the order of the mean free path or greater, the different centers 
(ri,r2 = ri ± (Tij) of the colliding pair in Eq. (|3.9p can be neglected (ri « r2) since the diameters of the particles 
are small compared to the mean free path at low density. As will be shown below, a nonzero distance between the 
particle centers gives rise to the collisional contributions to the transport coefficients, which are not present in dilute 
systems. These two modifications to fij result in the usual Boltzmann description for a granular mixture. The results 
obtained here therefore encompass earlier work on granular mixtures at low density [38*1 . In the elastic limit, aij — *■ 1, 
these equations become the Enskog theory for mixtures of dense molecular gases studied in Ref. [4§|. 

As happens for elastic collisions, the inelastic Enskog equation provides a semiquantitative description of the 
hard sphere system that neglects the velocity correlations between the particles that are about to collide (molecular 
chaos assumption). The Enskog approximation is expected to be valid for short times since as the system evolves 
corrections to the Enskog equation due to multiparticle collisions, including recoUision events ("ring" collisions) should 
be incorporated. The latter are expected to be stronger for fluids with inelastic collisions where the colliding pairs 
tend to be more focused. Therefore, some deviations from molecular chaos have been observed in molecular dynamics 
(MD) simulations [53, "ss", '56| of granular fluids as the density increases. Although the existence of these correlations 
restricts the range of validity of the Enskog equation, there is substantial evidence in the literature for the validity 
of the Enskog theory at moderate densities and higher restitution coefficients especially at the level of macroscopic 
properties (such as transport coefficients). In the case of molecular dynamics (MD) simulations, the Enskog theory 
compares quite well with simulations for the radial distribution function 53], the self-diffusion coefficient [57, 58], the 
kinetic temperatures of a binary mixture in homogeneous cooling state [59] , and the rheological properties of a mixture 
under simple shear flow (goI . [6]| . The agreement between MD and Enskog equation is good for moderate densities 
(solid volume fraction up to 0.15) and even conditions of strong dissipation (restitution coefficients > 0.7). For 
higher densities the a range is more limited but the Enskog theory still captures the relevant qualitative features. 
The Enskog transport coefficients for a monocornponent gas ^1] have also been tested against NMR experiments 
of a system of mustard seeds vibrated vertically [H, . The average value of the coefficient of restitution of the 
grains used in this experiment is a = 0.87, which lies outside of the quasielastic limit (a w 0.99). Comparison 
between theory and experiments shows that the Enskog kinetic theory successfully models the density and granular 
temperature profiles awa y fr om the vibrating container bottom and quantitatively explains the temperature inversion 
observed in experiments [63] . All these results clearly show the applicability of the Enskog theory for densities outside 
the Boltzmann limit and values of dissipation beyond the quasielastic limit. In this context, one can conclude that 
the Enskog equation provides a unique basis for the description of dynamics across a wide range of densities, length 
scales, and degrees of dissipation. No other theory with such generality exists. 



IV. MACROSCOPIC BALANCE EQUATIONS 



In the previous section, the Enskog assumption (|3.6p was used to obtain a closed set of kinetic equations (|3.8p for 
a moderately dense mixture of inelastic hard spheres. The result takes the form of nonlinear, integral-differential 
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equations for the distribution function /i, which contains information on a microscopic scale. In this section, this 
theory will be used to obtain the corresponding description on the macroscopic (or hydrodynamic) scale. First, 
the relevant macroscopic variables will be identified and defined. Next, the corresponding balance equations will be 
derived. Finally, expressions for the equations of state (pressure and cooling rate) and fluxes will be presented as 
integral expressions containing fi. 

The variables of interest for a macroscopic description of the mixture are the number densities for all species, 
{ui (r, t)} (or equivalently, the mass densities {pi (r, t) = rniUi (r, t)}), the total energy density, e (r, t), and the total 
momentum, p (r, t) . These are expected to be the s + 1 + d slow variables that dominate the dynamics for long 
times through a closed autonomous set of equations, the hydrodynamic equations. The reasoning behind this is 
that these are the densities for global conserved quantities in molecular fluids, and therefore have decay times set by 
the wavelength of the excitations. Long wavelength (space scales large compared to the mean free path) phenomena 
therefore persist at long times (compared to a mean free time) after which the complex transient microscopic dynamics 
has become negligible. For granular fluids, the energy is not conserved but is characterized by a cooling rate at long 
wavelengths. Still, this cooling rate may be slow compared to the transient dynamics and thus the energy remains a 
relevant slow variable. This is confirmed by MD simulations showing a rapid approach to this cooling law after only 
a few collisions [59] . 

These macroscopic variables will be referred to collectively as the hydrodynamic fields. They are defined without 
approximation in terms of moments of the distribution functions 

rii (r, i) EE y dv/j(r,v;i), i=l,..s, (4.1) 
e{r,t)=J2 J dv^m.v^Mr, v; t) (4.2) 

2 — 1 

p{r,t) = Y, f dvm,v/,(r, v; t) (4.3) 

The time dependence occurs entirely through the distribution function and hence is determined from the Enskog kinetic 
equations (j3.8p . However, rather than solving the kinetic equation to determine this complete time dependence it is 
useful for the purposes of deriving the simpler hydrodynamic description to first obtain the balance equations. These 
equations express the time derivative of the hydrodynamic fields in terms of local fluxes and sources due to collisions 
or the external force. These equations and the identification of the fluxes follow in detail from the form of the collision 
operators in ()3.2p as shown in Appendix |B] (in fact they are obtained there exactly from the first hierarchy equation 
(|3.ip without the Enskog approximation (|3.6p and hence are exact). The results for the balance equations are 

dtn, (r, <) + m~i V • j, (r, t) = 0, (4.4) 

S 

dte (r, i) + V • s (r, t) = -w (r, t) + ^ m-^F, (r) • j, (r, t) , (4.5) 

i=l 

s 

dtpp (r, t) + dr^t^p (r, t) = Y,n^ (r, t) (r) . (4.6) 

1=1 

The explicit expressions for ji, s, w and t^p are contained in Appendix iBl and not shown here since they are cast in a 
more convenient form below. 

The mass fluxes {ji (r, t)} , energy flux s (r, i), and momentum flux tp^ (r, t) describe the rate of transport of the 
hydrodynamic flelds through a given cross sectional area. They consist of parts due to pure convection and parts due 
to collision. To identify the convective (kinetic) parts, the local flow field U (r, t) is defined in terms of the momentum 
density by 

k 

p(r,0 = p(r,t)U(r,t), p (r, t) = ^ m,n, (r, i) , (4.7) 
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where the second equation defines the mass density. Also, the energy density is written in terms of the internal energy 
density Cq (r, t) in the local rest frame, plus the energy due to flow 

e (r, t) = eo (r, t) + ip (r, t) (r, t) . (4.8) 

In terms of U (r, t) the fluxes become 

j, (r, t) = p, (r, t) U (r, t) + jo. (r, t) , (4.9) 

S0 (r, t) - ( eo (r, t) + (r, t) (r, t) ) Uf3 (r, t) + Pp^ (r, t) (r, t) + qp (r, t) , (4.10) 

tp^ (r, t) = p (r, t) Up (r, t) t/^ (r, t) + Pp^ (r, t) . (4.11) 

The first terms on the right sides describe convective transport, while the diffusion fluxes joi (r,f) , heat flux q(r, t), 
and pressure tensor Pp^ (ri , t) describe the residual transport for each fluid element in its local rest frame. Before 
giving their forms more explicitly, it is instructive to insert (|4.9|l - (|4.1ip into ( I4.4p - (|4.6p to get the equivalent form 
for the balance equations 

An^ + WiV-U + m^V-joi = 0, (4.12) 

S 

Dteo + {eoS^p + P^p) d^^Up + V ■ -w (r, t) + ^ mr^F, (r) • jo, (r, t) , (4.13) 

i=l 

s 

pDtUp + dr^P^p = X! ' (4-14) 

where = 9t + U • V is the material derivative. 

The independent hydrodynamic fields are now {ni (r, i)}, cq (r, i)i U(r, i). The remaining quantities in the 
balance equations are the energy loss rate w (r,t), the mass fluxes {joi (r,i)}, the heat flux q(r,t), and the pressure 
tensor Pp^ {vi^t). These quantities, which are defined in terms of the distribution functions, are obtained by the 
explicit forms for j^, s, w, and t^p given in Appendix IB] together with Eqs. (|4.9p - (|4.1ip . 

Specifically, the energy loss rate is due to inelastic collisions 



\3 , 



xe{d- ■ gi2)(5^ • gi2)'^/jj(ri,vi,ri +cr,j,v2;t), (4.15) 
whereas the diffusion flux arises from convective (kinetic) transport 

joi(ri,t) = mj y dviVi/j(ri, vi;t), (4.16) 

where Vi = vi — U(r,t) is the velocity in the local rest frame. The heat flux has both "kinetic" and "coUisional" 
transfer contributions 

q(ri,i)EEq'=(ri,i)+q=(ri,t), (4.17) 

with 

q'=(ri,t) = ^ J dviim,Fi2Vi/,(ri,Vi;i), (4.18) 



8 



k 

q"" (ri, t) = X! ^ + a^j)mjfi,jafj J dvi J J da e{a ■ gia) 



^2 



^{S--gi2) [(1 - a^j) {^iji - ^lij) {^ • gi2) + 4ct ■ G.y) 

XCT / (ix/ij (ri - Wy , vi,ri + (1 - x) (T,y , V2;t), 
Jo 



(4.19) 



where Gy = /ijjVi + /ijiV2 is the center-of-mass velocity. 

Similarly, the pressure tensor has both kinetic and coUisional contributions 



where 



P^p (n, i) = (n, + P^p (n, , (4.20) 



Pl;^ (ri,i) = dvim,VipVi-J,in,v,;t), (4.21) 



4=1 



-P7/3(l"l,i) = ^ ^ mjfJ-tA'^ + J C^Vl y dv2 /"dS^e(ff • gi2)(? • gl2)^ 

X(T/3(T7 / da;/ij(ri - 2:<Tij,vi,ri + (1 - x)<Tij,V2;t). 
Jo 



(4.22) 



Equations (|4.12[) - (|4.14p together with the definitions (|4.15p ~ (|4.22|) represent the macroscopic balance equations for a 
granular mixture, without restrictions on the densities or degrees of dissipation. In the case of a three-dimensional 
system (d = 3), the above equations reduce to previous results ^3] derived for hard spheres. When the approximate 
form (|3.6p is used in the first hierarchy equation and in these expressions for the cooling rate and fluxes, the Enskog 
theory results. 

For historical consistency with the usual constitutive equations for a ordinary fluid, the temperature T (r, t) is used 
in the following instead of the internal energy density eg (r, t), with the definition 

eo(r,i) = ^n(r,i)r(r,t). (4.23) 

As a definition, this amounts only to a change of variables and there are no thermodynamic implications involved in 
the use of this temperature for a granular fluid. The corresponding hydrodynamic equation for T (r, t) follows directly 
from (|4l3l) 

d d * * 

-n {Dt + 0T + P^pdr^ f/^ + V • q - -T ^ m" V • jo, = ^ m-^F, • jo,. (4.24) 

i=l i=l 

To obtain these results the continuity equation has been used 

Dtp + pV -V = 0. (4.25) 

This follows from the definitions of p and U and the conservation laws for the {rii (r, i)} . A related consequence is 



I]jo. = 0, (4.26) 



so that only s — 1 dissipative mass fluxes are independent. Finally, the "cooling rate" C, has been introduced in (|4.24p 
by the definition 



dnT 2dnT 



i— ^ (l - afj) m,p,ji(jf^ M dvi / dv2 / rfS, 



xe(? • gi2)(5^ • gi2)Vu(ri,Vi,ri + (Ty,V2;t). (4.27) 
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V. CONCEPT OF A NORMAL SOLUTION AND HYDRODYNAMICS 

The form of the equations of state and fluxes given in the previous section, (|4.16|) - (|4.22p and (|4.27[) . are cast as 
functionals of the distributions {/i}, which depend expUcitly on space and time. As a result, the macroscopic balance 
equations are not entirely expressed in terms of the hydrodynamic fields, and thus do not comprise a closed set of 
equations. If these distributions can instead be expressed as functionals of the hydrodynamic fields (norma? solution), 
then C (r, t), {jpi (r, t)}, q (r, t), and P^^ {^i,t) also will become functionals of the hydrodynamic fields through (|4.16p - 
(|4.22p and (|4.27p . Such expressions are called "constitutive relations". They provide the missing link between the 
balance equations and a closed set of equations for the hydrodynamic fields alone. Such a closed set of equations 
defines "hydrodynamics" in its most general sense. 

It is seen, therefore, that any derivation of hydrodynamics proceeds first by construction of normal solutions to the 
kinetic equations. More precisely, a normal solution is one whose space and time dependence occurs entirely through 
the hydrodynamic fields, denoted 

/,(ri,vi;t) = /,(vi |{y;3(ri,i)}), (5.1) 
where {y^ (ri,t)} denotes generically the set of hydrodynamic fields 

y0^{T,U,{n,(r,t)}}. (5.2) 
Therefore, the space and time derivatives of the kinetic equation are given by 

{dt + vi . Vr) /.(vi I {yf, (n, i)}) = / I W)}) + . Vr) (r; t) . (5.3) 

J Syr, (r; t) 

Furthermore, the balance equations for the hydrodynamic fields (|4.f 2p - (|4.f 4p can be used to express dtyrj (r; t) in 
(|5.3p in terms of space derivatives of the hydrodynamic fields. For such a solution for fi, Eqs. (|4.16p - (|4.22p and (|4.27p 
give directly by integration the desired constitutive relations. 

The determination of /i(vi | {yp (ri,t)}) from the kinetic equations (|3.8p is a very difficult task in general, and 
further restriction on the class of problems considered is required at this point to make progress. Any functional of 
the fields can be represented equivalently as a local function of the fields and all of their gradients. In many cases, 
gradients of high degree are small and may be negligible so that the normal distribution becomes 

/,(vi I {y/3(ri,i)})^/,(vi;{2/^(ri,0,Vriy/3(ri,t),---}) (5.4) 

This representation does not imply that the low degree gradients are small, and fi may be a non-linear function of 
the relevant gradients. This occurs in many important applications for granular fluids 40]. In the limiting case where 
the low-degree gradients can be controlled by boundary or initial conditions and made small, a further Taylor series 
expansion can be given 

/.(vi I {y0 (ri,0}) - /^(vi; {y^ (ri,i)}) + /P^(vi; {yp (ri,t)}) + • • • 

/f^(vi; {y/3 (ri,i)}) -I- Y,,(vi; {y^ (ri,i)}) • Vr^j/a (ri,0 • • • 

(5.5) 

It follows that the leading order distributions have the exact properties 

(r, t) = J dwfl''^ (v; {yp (r, i)}), * = 1, -s, (5.6) 

^n{v,t)T{v,t)^j2 j dvim,yVi°^(v;fe(r,i)}), (5.7) 
p (r, t) U (r, t) EE ^ y dvm,v/f ) (v; {yp (r, t)}), (5.8) 

i—l 

and the corresponding moments of all higher order terms in (j5.5p must vanish. Generalization of this type of gradient 
expansion for the normal solution to include a class of nonlinear gradients in the reference state has been discussed 
recently [ssLIbgI. 
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As is standard for molecular gases, the gradient expansion will be taken with respect to the reference local HCS, i.e. 
that resulting from the neglect of all gradients in the functional but evaluated at the value of the fields at the chosen 
point and time {?/^ (ri, t)}. This point is crucial in our analysis since most of the previous results have taken elastic 
Maxwell distributions as the base state. Note that in the CE method the form of /j^"-* comes from the solution to 
the kinetic equation to zeroth order in gradients and cannot be chosen a priori. Accordingly, fl^\'Vi', {yp (ri,t)}) — > 
f-^\Vi]{yf3{ri,t)}), where F=|v — U(r,t)| is homogeneous and isotropic with respect to its velocity dependence. 
This symmetry implies that the leading (zero) order contributions to (|4.16p and (|4.17p for the vector fluxes {joi (r,t)} 
and q(r,t) must vanish, and this contribution to the pressure tensor P^p must be isotropic (proportional to S-y^). 
Similar symmetry considerations to the first order contribution (linear in the gradients) determines the exact structure 
of the constitutive equation to this order. Based on these symmetry considerations, the constitutive quantities are 
known to take the forms 

C (r, t) ^ C(°) {{yp (r, t)}) + Cu {{y0 (r, t)}) V • U (r, i) , (5.9) 

* n ■ (r t) 

joj(r,i) -^m^TTij ^ ' Dij{{yp{v,t)})V\nnj{Y,t) 

S 

-p (r, t) Dj {{yp (r, t)}) V InT (r, i) - ^ Df^ {{y^ (r, t)}) F, (r) , (5.10) 
q(r,i) -A({y^(r,t)})VT(r,t) 



^ (T^ (r, t) {{yp (r, t)}) V Inn, (r, t) + {{yp (r, t)}) F, (r)) 



(5.11) 



P^A (r, t) ^ p {{yp (r, t)}) - V {{yp (r, t)}) (^dr^Ux (r, t) + dr, (r, t) - • U (r, t) 

-n{{yp{v,t)})V -Viv.t). (5.12) 

The unknown quantities in these constitutive equations (|5.9p - (|5.12p include the cooling rate C^"-' ({y/j (r, t)}), the 
hydrostatic granular pressure p[{yp{v,t)}), and the transport coefScients Qu, Dij, Df , Dfj, Dq,ij, Lij, 77, and k. 

These quantities can be expressed as explicit functions of the hydrodynamic variables once /j^^'' and f^^"^ are known. 

The equations governing the solution of f^^'^ and f^^^ are found using the CE method, as described below. 

VI. CHAPMAN-ENSKOG NORMAL SOLUTION 

The CE method is a procedure for constructing an approximate normal solution. It is perturbative, using the 
spatial gradients as the small expansion parameter. More precisely, the small parameter is Knudsen number (Kn), 
defined as the gradient of the hydrodynamic fields relative to their local value times the mean free path. This means 
that the conditions for the solution are restricted to small variations of the hydrodynamic fields over distances of the 
order of the mean free path. In the presence of an external force it is necessary to characterize the magnitude of this 
force relative to the gradients as well. Here, it is assumed that the magnitude of the force is first order in perturbation 
expansion. This allows comparison with the results of Ref. [i^ for the elastic case. 

The perturbation is carried out by considering the Enskog kinetic equations successively at each order in the 
gradients. As described below, the zeroth order equation is first obtained for f^^\ Next, the first order equation for 
/f^"* is obtained. This expansion leads to integral-differential equations for the determination of /f"'' and /P"*, which 
are solved explicitly in the follow-on paper [50|. 

As detailed in Appendix [Cl to zeroth order in the gradients, the kinetic equation (|3.8p becomes 

s 

{d,T)dTf^'\vi;{yp{n,t)}) h I ' (6-1) 
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where 



t(0) 



dv2 / daQ{a ■ gi2)(CT • gi2) 



iVi'; {yp (ri, ) (F^"; {y^ (n, t)}) 



/f)(Fi;{y^(ri,t)})/)"^(T/2;{y^(ri,t)}) 



(0), 



(6.2) 



All spatial gradients arc neglected at this lowest order. Equation (|6.ip determines the velocity dependence of 

fl^\Vi] {yp (ri, t)}); the space and time dependence is local and entirely through the fields {yp {vi,t)} at the space 

and time point of interest. This has been exploited by writing the time dependence of /I"' in terms of the time 
dependence of the fields, and recognizing that all time derivatives of the latter are proportional to space gradients, 
except the temperature, through the balance equations 



dtn, = 0, dtT - -(("^T, dtV 
Here, C^"^ is the cooling rate ()4.27|) to zeroth order in the gradients 



0. 



(6.3) 



1 



2d7iT 



^ (l - aj^) m,^j,x[f i<J^j■,{n^ (ri, i)}) afj ^ J d^i J dvs J da 



xe(? • g^2){^ ■ si2)'f^'\Vu{y0 (ri,t)})/r(y2; {yp in^t)}) 



(0)/ 



(6.4) 



)(' 

Similarly, the functional dependence of Xij^ (ri, ^2 \ {^i}) on the compositions to zeroth order in the gradients has the 
functional dependence on the densities replaced by {rii} —f {ui (ri, i)}, at the point of interest. The result is transla- 



(0) 
At J 



|ri — r2| ; {rii {ri,t)}) , a function of the densities. Finally, 



tional and rotational invariant x^j^ (ri,r2 | {rii}) 

gradients in the distribution functions of the collision operators must be neglected, e.g. /^"^(vi; {yp (ri + cr^ , t)}) 
/(°)(vi;{y/3(ri,t)}). 

A further simplification of these equations for the lowest order distribution functions occurs when they are written 
in terms of the corresponding dimensionless forms {(t>i} 



friV;{y,})^n,v~'^{T)cl>,{V*;{n*}), 



(6.5) 



with the definitions 



V* 



voiT)' 



V m 



1 ^ 



(6.6) 



The solution depends on the flow field only through the relative velocity V. Furthermore, since there is no external 
energy scale the temperature can occur only through the scaling of the dimensionless velocity through the thermal 
velocity vq{T). Equation (|6.1[) now takes the dimensionless form 



(6.7) 



where Vv* = d/dY* and 



'■(0) 



AO)* 



■''ij ' 



s 1 ^ 



(6.8) 



i=l 



i=l 



The solution to this equation is a universal function of the magnitude of the velocity V* and is otherwise independent 
of the temperature and flow field. For a one component fluid it is independent of the density as well. However, for 
mixtures it is parameterized by the dimensionless species densities through the factors Xij^ {^ij] {ri*}) ■ Equation (|6.1[) 
has the same form as the corresponding dimensionless Enskog equations for a strictly homogeneous state. The latter 
is called the HCS. Here, however, the state is not homogeneous because of the requirements (I5.6p - (l5.8p . Instead it is 
a local HCS. As said before, an important point to recognize is that the occurrence of this local HCS as the reference 
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state is not an assumption of the CE expansion but rather a consequence of the kinetic equations to lowest order in 
the gradient expansion. 

The analysis to first order in the gradients is similar and the details are given in Appendix [C] The result has the 
form (EI 



f(i) 



A (V) • V InT + ^ Bl (V) • V In) 



+Ci,7r, (V) - ( + dr^U^ - -(5^„V • U 



-I?.(V)V.U + ^£nV)-F, 



(6.9) 



The contributions from the flow field gradients have been separated into independent traceless and diagonal com- 
ponents, as follows from fluid symmetry. The velocity dependence of the gradient contributions is contained in the 
functions A (V) , (V) ,Ci^^^ (V) ,1?^ (V), and £\{y). The kinetic equations determine these functions as the 
solutions to the integral equations 



(6.10) 



(6.11) 



(6.12) 



(6.13) 



(6.14) 



The linear operator C is given by 



(£X), = ic(°)Vv(VXO + (XX),, 



(6.15) 



AO) 

ij 



VI I X,,/ 



(0) 



7(0) 



(6.16) 



and the inhomogeneous terms are defined by 



Ar.-r (V) = ^V.Vv • (v/f 



.7=1 



(6.17) 



(0)n 



f(o) 



(6.18) 
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C^,,p (V) = \ (v.dv.ff'^ + Vpdvjl'^ - -^5p,Y ■ Vv/f ^) 

+ 1 E fe^7K/i°^] + J^vAdvjf] - j,6p,JC,,AdvJ^°^]) , (6.19) 



A(V) = iv-Vv/f -i(Ct/ 



nTd 



f(0) 



(0) 



ir 



(0) 



(6.20) 



E;,V,.-(Vv/;»>)i-(..,-!^), (6.1) 

where the operator JCij^^[X] is defined by Eq. (jC16p . 

This completes the construction of the normal solution to the revised Enskog equations up through first order 
in the gradients. Equation ()6.7p determines the /I"' through the definition (16. 5^ : solution to the linear integral 
equations (|6.10p - (|6.14p determines the f^^'^ through the definition (|6.9p . The unknown fiuxes and cooling rate of the 
hydrodynamic equations can then be calculated with these solutions. This is made explicit in the next section. 



VII. CONSTITUTIVE EQUATIONS AND TRANSPORT COEFFICIENTS 

The forms for the constitutive equations to first order in the gradients are given by Eqs. (|5.9p - (|5.12p . The explicit 
representations for the coefficients in these equations are given in terms of the solutions to the integral equations for 
ff''^ and /j^^^ of the previous section. Details of the simplification of these expressions in terms of ff''^ and f^^^ are 
given in Appendix iFl and only the final results are presented here. 

Recall that the results below are based on the assumption that the external force is of the same magnitude as /f^"*; 
a force of different magnitude would result in different constitutive relations. 



A. Cooling Rate 

The cooling rate is calculated from Eq. (|4.27p . resulting in the form (|5.9p 

C - C^°^ + Cc/V • U, (7.1) 

with 

<<»' - ^ i: (1 - 4) / / .v./.'«'(v'.,/f (7.2, 

1,3 = 1 ■> 

and 

(i+2 , 2\ (0) d T,(o) 

+lf t (1 - ^>^S<' I / <iv, .?,/™(V.)I>,(V,,. (7.3) 

l,j = l 

The constant Bn is defined by 
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Also in (|7.3p the species temperatures have been defined by 



(7.5) 



There is no special significance to these quantities other than naming the integral on the right side, which is a specified 
function of the hydrodynamic fields {rii} and the global temperature T through fj^^\ 



B. Mass Fluxes 

The mass fluxes are determined from the definition of (|4.16|) leading to the form (|5.10p to first order in the gradients 

s 

jo, ^ - ^ m^mJ^D,J\/ In - pDiV InT - } ] D^^F, . (7.6) 
The transport coefficients are identified as 



_ rmmj^D.jV In rij ~ pDfV In T - ^ D[^Fj 



D 



A, = 



pd 



nijnjd 



J dvY ■ B{ (V) 



A^ = -^ / rfvV.^Hv), 



(7.7) 
(7.8) 
(7.9) 



C. Energy Flux 

The energy flux to first order in the gradients is given by (|5.1ip 

s 



There are both kinetic and coUisional transfer contributions according to Eq. (I4.17p . q 
contributions to the transport coefficients are identified as 



(7.10) 



q + q"^. The kinetic 



- E - E / rfvim,y^v . A (V) , 

4=-^ / rfvim,^^V.£^(V). 



(7.11) 



(7.12) 



(7.13) 



For convenience below, the partial thermal conductivities have been introduced in Eq. (|7.1ip . The collision transfer 
contributions are obtained from (|4.19p to first order in the gradients. These are calculated in Appendix iFl with the 
results 



8S2 



2 + d 



^A^-(d + 2) 



T, 



(0) 



rriimjT 



{2pij - Pji)pD 



2 T^°^ 
— A,^ + (d + 2)^— 



iriT 



pDj 



(7.14) 
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" 1 



{d + 2] 



prrii 



1 J J T~i 

^^2 P 

1 ^ rriip 



2r'T 



(7.15) 



^Ij = ^zi^+^w)'^P^^^p(^^pXfJ {^Bi{l-a^p){^l^p- ppi) 



p=l 



— 11^ + (d + 2) 



t: 



(0) 



-D 



p] 



8B2 

? 

d + 2 



'^t^pi J- k 



P3 



(d + 2) {2fiip - ppi) 



(0) 



miirLp 



P3 



(7.16) 



where the coefficients Cj^ and Cf^^ are given by Eqs. (|F27p and (|F28|) . These expressions also depend on the transport 
coefficients of the mass fluxes, Dj , Dij, and Dfj given by Eqs. (|7.7p . (|7.8p . and (|7.9p . respectively, and on the kinetic 



contributions Af , , and L^j . 



D. Momentum Flux 

The pressure tensor is evaluated from Eqs. (|4.20p - (|4.22|) . To zeroth order in the gradients, one gets the pressure p 



({n J , T) = ^ / ({n,} ,T)+p^ ({n J , T) = ip(f + ip(f , 



where 



(7.17) 



(7.18) 



(0) 



Similarly the shear viscosity is 77 = rj'' + rj'^ where 



i=l 



id + 2){d 



(7.19) 



(7.20) 



2Bi ^ . .(o)^_d^fc , 



(c^ + 2) 

Finally, the bulk viscosity is k = k'^ + k'^ where 

S3 (rf + 1) 



d + 2 



(7.21) 



= 0, k'^ 



2d2 



^ m,^., (1 + a.,) xgV^+' / dvi | dv^/f ^ (Vi)/f ^ (V2)5i2. (7.22) 
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VIII. DISCUSSION 

The most complete and accurate description of mixtures for ordinary fluids is based on the revised Enskog kinetic 
equations for hard spheres. The expUcit construction of solutions to those equations by the CE expansion to first 
order in the gradients was given more that twenty years ago in Ref. [i^. These solutions, together with the macro- 
scopic balance equations obtained from the kinetic equations, provide a self-consistent derivation of Navier-Stokes 
hydrodynamics for mixtures and the identification of expressions for all the Navier-Stokes parameters (equations of 
state, transport coefficients). In the context of the chosen kinetic theory, the analysis and the expressions for these 
parameters are exact. At this formal level questions of principle could be addressed, prior to the introduction of 
subsequent approximations for practical evaluations. For example, it was shown that application of the analysis to 
the original and revised Enskog theories leads to qualitatively different Navier-Stokes hydrodynamics, only one of 
which is consistent with irreversible thermodynamics. Since no approximations were involved this was sufficient to 
reject the Enskog kinetic theory in favor of its revised version [isj . 

The present work is simply an extension of that in Ref. [ioj to inelastic hard sphere granular mixtures. Modification 
of the collisions to account for inelasticity leads to significant differences from ordinary fluids in detail, but the formal 
structure of the CE expansion remains the same. Similarly, granular Navier-Stokes hydrodynamics results exactly 
from the CE solution to first order in the gradients and the corresponding modified balance equations. The form of 
these hydrodynamic equations and expressions for the transport coefficients are exact, as in the ordinary fluid case. 
The primary motivation for this analysis is to provide the basis for practical applications, as noted in the Introduction, 
and described in the following paper. However, at the formal level, important fundamental questions can be addressed 
and clarifled as well. 

The existence of hydrodynamics for granular fluids has been questioned, due to the many known differences from 
ordinary fluids: there is no equilibrium or even stationary reference state; the temperature is not a hydrodynamic 
field (failure of energy conservation) , or conversely, multiple temperature fields could be required for mixtures (failure 
of equilibrium state equipartition for the corresponding granular HCS). In the end, qualitative discussions must be 
resolved by controlled analysis. Here, the validity of the RET for some range of densities and degree of dissipation has 
been assumed as a mesoscopic basis for possible macroscopic dynamics in a granular mixture. As shown in the text, 
sufficient conditions are the macroscopic balance equations (verified) and a normal solution to the kinetic equations. 
The normal solution is defined in terms of a chosen set of hydrodynamic fields, and the question of hydrodynamics 
reduces to its existence. The details of the Appendices give the explicit construction of this solution to first order 
in the gradients, together with a proof of the existence of solutions to the associated integral equations. It can be 
concluded from this that a closed set of hydrodynamic equations for the species densities, flow velocity, and a single 
temperature exist for sufficiently small gradients. 

This conclusion is consistent with the observations that the reference state is not equilibrium, depends on the cooling 
temperature, energy loss can be large at strong dissipation, and the kinetic temperatures of species are different. None 
of these facts compromises implementation of the CE expansion for solution to the kinetic equation. The parameters of 
the resulting Navier-Stokes equations incorporate such effects through the integral equations that determine them, and 
their dependence on the time dependent fields. This in turn affects the solutions to the Navier-Stokes equations under 
different physical conditions, and is responsible for some of the observed peculiarities of granular fluids. Clearly, 
it is important to get the details of the Navier-Stokes equations accurately before concluding that any observed 
experimental phenomenon is hydrodynamic or not. This is another primary motivation for the present work. 

These details entail solution to the equation for the reference state and solution to the integral equations for the 
transport coefficients, to determine them as functions of the hydrodynamic fields (temperature, flow field, and species 
densities) and the system parameters (restitution coefficients, masses, particle sizes). There has been considerable 
study of the reference state, as an expansion about a Gaussian for relatively small velocities (asymptotic forms for 
large velocities are known as well). The integral equations can be solved approximately as truncated expansions in a 
complete set of polynomials with Gaussian weight factors. For ordinary fiuids the leading approximation is generally 
quite accurate, and the following paper gives its extension to the granular mixture. Still, there are open questions 
about this approximation for strong dissipation and large mechanical disparity (e.g., mass ratio). Previous results 
obtained for granular mixtures at low-density [63, and for the shear viscosity of a heated granular mixture at 
moderate density [63 | have shown the accuracy of the above approximation, even for strong dissipation. 

An accurate solution to the integral equations will predict the transport coefficients as functions of the dissipation. 
There is only one correct result for this dependence, given by the formulas obtained here. However, its measurement 
in a given experiment or simulation can entangle and affect this dependence of the transport coefficients due to higher 
order gradients beyond the Navier-Stokes limit. It may be tempting to compare experimental or simulation data to 
a corresponding Navier-Stokes solution, adjusting the transport coefficients for a best fit and reporting these as the 
"measured" values. This can be misleading for granular fiuids under conditions where the size of the gradients increase 
with the degree of dissipation. For such states, strong dissipation can require additional terms in the constitutive 
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equations beyond those of Navier-Stokes order [i^, [sl, [T^I ■ This does not mean that the results obtained here are not 
correct at strong dissipation, only that they must be distinguished carefully from other effects of the same order. A 
careful tabulation of the Navier-Stokes results given here (e.g., via Monte Carlo simulation) is required for an accurate 
analysis of experiments of current interest. It is an interesting new feature of granular fluids that hydrodynamic states 
beyond Navier-Stokes order may be the norm rather than the exception. 
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APPENDIX A: RET AND SPATIAL CORRELATIONS AT CONTACT 

In the case of ordinary fluids, the Enskog approximation can be understood as a short time, or Markovian approx- 
imation. This follows if the initial distribution has the Enskog form 

/y(ri,vi,ri - cr,y, V2;t = 0) = fi{ri,vi,t = 0)/j(ri - crij,V2;t = 0)xy (ri,ri - (T^j \ {rifc}) . (Al) 

In fact, this is a quite plausible class of initial conditions since correlations in that case are generally induced by the 
interparticle structure that is independent of the velocities. Then at finite times, it is assumed that /ij(ri,vi,ri — 
(Tij , V2 ; t) becomes a functional of fi 

/ij(ri,vi,ri - cTij,X2]t) = ^'y(ri,vi,ri - cry,v2;t | fi{t)). (A2) 

The Enskog approximation corresponds to evaluating this functional at t = 

/jj(ri,vi,ri - cry,v2;t) Fjj(ri,vi,ri - cr„-,v2;t = | f^{t)). (A3) 

Thus for the special class of initial conditions the Enskog approximation is asymptotically exact at short times, and 
assumes that the generator for dynamics at later times is the same as that initially. This idea provides a simple mean 
field theory for particles with continuous potentials of interaction, but is more realistic for hard spheres where there 
is instantaneous momentum transport at the initial time. The presence of inherent velocity correlations for granular 
fluids suggests that the form (|A1[) is less justified than in the ordinary fluid case. However, it is noted that velocity 
correlations are present for any nonequilibrium state even with elastic collisions and it is known that the Enskog 
equation still provides a good approximation in these latter cases. 

An important exact boundary condition for hard spheres is given by f5ji| 

0(S^ • gi2)/y(ri,vi,ri - (T.^,W2]t) = a~.^fe~-^e(-CT • gi2)/ij(ri,vi,ri - aij,\2]t). (A4) 

This equation implies that the distribution of particles that have collided is the same as those about to collide, but 
with their velocities changed according to the collision rules. In general the two particle distribution function can be 
written as 

/y(ri,vi,ri - (T,j,V2;t) = e(-5^ • gi2)/»j(ri, Vi, ri - cr,j, V2; i) 

+e(CT • gi2)/ij(ri,vi,ri - (Tij,v2;t). (A5) 

If the Enskog approximation (jA3[) is introduced in the first term on the right side of (jASp . then the corresponding 
approximation on the right side of Eq. (|A5p gives the approximate two particle distribution function at contact as 

/„ (ri,vi,ri - cr„-,V2;i) e(-5^ • gi2)Xjj (i"i, i"i - 0-^^ | {rife}) /i(ri, vi; t)/j (ri - 0-^^ , V2; t) 

+ a^r.'^b'^'^ei-^ ■ gi2)Xij iri,ri - cr,j \ {rifc}) /,(ri, vi; i) 

x/j(i"i - crji,V2;t) 
= e(-o- • gi2)x»j (ri,ri - \ {uk}) Mri,vi;t)fj{ri - (Tij,V2]t) 

+a~j'^e{a ■ gi2)Xij (ri,ri - cr.j \ {rifc}) /,(ri, v'/; t) 

xfjin- cr,j,v'^;t). (A6) 
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Since v" and are functions of both Vi and V2 there are velocity correlations on the complementary hemisphere, 
even when they are neglected on the precoUision hemisphere. 

An important consequence of (|A6P is the relationship of Xij (ri,ri — cTij \ {ni}) to the pair correlation function 
gij (ri,ri - cr,., | {rifc}) defined by 

Hi (ri) rij (ri ~ (Tij) gij {ri,ri- (Tij \ {nk}) ^ J '^^i y rfv2 /^^(ri, vi, ri - cr^ , V2; t) (A7) 
Use of the approximation (|A6p gives the result 

Hi {ri) rij {ri - (Tij) gij {ri,ri - (Tij \ {uk}) = 1 + «»3 (n, n - crij \ {uk}) ( dwi [ dv2 

ai] J J 

xe(-CT • gi2)/j(ri,vi;t)/,(ri - (T,j,W2]t), 

(AS) 

where a change of variables has been made in the integration of the second term in (|A6P 

j dvi j dv2X(v'/, v^') - a,j J dV; J dVi X(v'/, v^')- (A9) 

For a uniform system, gij (ri,r2 | {nk}) gij (|ri — i"2| ; {nk}) and this expression reduces to 

1 + a ■ ■ 

g^3 {<Jtj;{nk}) = — — -Xij {^fc}) (AlO) 
zaij 

Equation (jAlOp is the result quoted in the text and provides the interpretation for Xij (ri,ri ~ crij I ■ For 

elastic collisions Xij {'^ijl {iT-k}) is indeed the pair correlation function at contact. The Enskog theory in that case 
takes Xij (ri,ri — | {uk}) to be the pair correlation function for an equilibrium nonuniform fluid whose densities 
are equal to those for the actual nonequilibrium state being considered. This assumption is based on the fact 
that structural correlations for hard spheres are entirely due to excluded volume effects which should be similar for 
equilibrium and nonequilibrium states. It is reasonable to extend this choice for Xij (ri,ri — cTy | {uk}) to granular 
fluids as well. Its accuracy can be judged by measuring (via MD simulation) the pair correlation given by (|A10[) with 
this choice on the right side. This has been done for the one component fluid, indicating reasonable results over a 
range of values for the restitution coefficient . 

APPENDIX B: BALANCE EQUATIONS AND FLUXES 

The macroscopic balance equations follow from the definitions (|4.ip - (|4.3p and the first hierarchy equation p.ip 

5tn, + Vri -ydvivi/, = y"dviC„ (Bl) 

dte + Wr, / dvi^TUivlvif, - ^F^ ■ I dvivi/j = ^ / dvi^m^vfQ, (B2) 

i=l •' i=l •' i=l •' 

dtPp + dr^^ ^ j d-viTJirViyViiifi - ^ UiFifj J dvirriiVifsCt. (B3) 

i—l i—1 i—1 

The integrals over the coUisional contribution C'i are analyzed below with the results 

dviC, = 0, (B4) 

J dvi^niiv'^Ci = -V ■ - w, (B5) 



i=l 
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/ d-vinriiVipCi = -dr^t^p. (B6) 

i=l 

Use of these expressions in ()Bip - (|B3|) gives the balance equations (I4.4p ~ ()4.6|) of the text with 

ji = rrii / dvivi/j, (B7) 



s 



dvi^m,vfvif, + s^, (B8) 



i—1 

The terms s"^, and f^^ arising from the coUisional contribution d are identified by further analysis of the left 
sides of (|B5|) and (jB6p . To do so consider the general expression for some arbitrary function ipi (vi) 

J dw^ij,Ci = ^^tf^ j dvi J d^2 J d5^e(a-gi2)(?-gi2)V'i(vi) 



-2 



^ [%-^/y(ri,v'/,ri - (Tij,W2]t) - /y(ri,vi,ri +cry,V2;t)] . (BIO) 

The restituting velocities are functions of the given velocities, v" — v"(vi,V2), = V2(vi,V2), defined by p.3p . 
These relations can be inverted to get 

vi = v" - /iji (1 + ay) (ct • gi2)CT, V2 = V2 + (1 + ajj) (ct • gia)^. (Bll) 

Therefore, in the first term of (|B10p it is possible to change integration variables from (ividv2 to dw'ldv2, with a 
Jacobian aij to get 

y dvi j dw2 j daQ{a ■ gi2)(CT • gi2)'(/'i (vi) v", ri - cT,j,V^;t) 

= y dv'/y dv'^l d5^e(-5^-g'/2)(-a-g'/2)^.(viK,vJ,'))/.j(ri,v'/,ri-<T,„vi';t) 

= y"rfvi ^£^^2 y"dCTe(CT • gi2)(CT • gi2)'(/'i (v'l (vi, V2)) /ij(ri, vi,ri + <Tij, V2;t), (B12) 

where use has been made of • gi2) = —aij{a- ■ g'/2)- In the last line the dummy variables (v'/,V2) have been 
relabelled (vi,V2), and a change of integration from a to —a has been performed. Accordingly vi (v'/jVj) has been 
relabelled v[ (vi,V2) with (jBll|) becoming in this notation 

v'l = vi - fiji (1 + aij) {a ■ gi2)5^, V2 = V2 + (1 + aij) {a ■ gi2)o^. (B13) 



This is the direct scattering law, which differs from the restituting scattering law (|3.3|) for 7^ 1. With this 
transformation the integral (|B10p is 

/ dviV'^a = E^y^'/ '^^i / / d^e(a.gi2)(^ -gia) 

X (V'i(v'i)-'0»(vi))/zj(ri,vi,ri+o-,j,v2;i). (B14) 

The special choice -0^ (vi) = 1 proves (jB4p above. 
Next, consider the sum of (|B14p over all species 

^ j dviV'^C, = ^ y dvi j dv2 y d^e(?-gi2)(S-gl2) 

X [V'j (v'l) - -0J (vi)] /^^^(ri,vi,ri +(Ty-,V2;t) 

" ^ XI "^y^^/ '^'^i / ^^"^2 / rfse(5^-gi2)(5^-gi2) 

X {[V'i (v'l) - (vi)]/ij(ri, vi,ri + (Ty,V2;i) 

+ [V'j (V2) - '(/'i (V2)] /ji(ri, V2,ri - cr.y, vi;t)}. (B15) 
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The second equality is obtained from the first by taking half the sum of the first plus an equivalent form obtained by 
interchanging vi and V2, interchanging i and j, and changing 5 to —a. To simplify this further, note the relation 
/3i(ri, V2,ri - cTy , vi;t) = fij (ri - cTj^, vi, ri, V2; i) and arrange terms as 

i=l jj = l 

X {[i^i (v'l) + ipj (V2) - V'i (vi) - ipj (V2)] /y(ri, vi, ri + cr^ , V2; i) 

+ bl^i (v'l) - ipi (vi)] [/y(ri,vi,ri + (Tij,V2;t) - fijiji - cry , vi, ri, V2; t)]} . 

(B16) 

The first term of the integrand represents a coUisional effect due to scattering with a change in the velocities. The 
second term is a collisional effect due to the spatial difference of the colliding pair. This second effect is called 
"coUisional transfer" . It can be written as a divergence through the identity 

/ij(ri,vi,ri +<Tij,V2;t) - /y (r 1 - cTy , vi, ri, V2; 

d 

= J dx— /jj(ri - xcTjj, vi, ri + (1 - a;) (Ty , V2;t) 

= Vri-o-jj- / dx fij{ri- xcrij,wi,ri + {1- x)cTij,V2]t). 
Jo 



(B17) 



Using the identity (|B17p . Eq. (jB16p can be finally written as 

J d^i^p^a = ^ E 4^'/ / J rf^e(s-gi2)(^-gi2) 



X {[^^ (v'l) + V-j (V2) - ipi (vi) - (V2)] /y(ri,vi,ri + (7,^,^2; t) 



+ 



Vri • (Tij [ipi (v'l) - Ipi (Vi)] J dx fijiji - XCTij, Vi,ri + {1-X) (Tij,V2;t)^ . 



(B18) 



Now, apply this result to the case ipi = rriiVi. Since the total momentum is conserved in all pair collisions, 
tpi (v'l) + V'i (V2) - i^i (vi) - tpj (V2) = for this case and (|B18P gives (jB6| with 

= / / '^'^2 / dCTe(CT • gi2)(ff • gl2)CTT, 

X (miv[i^ - miVi(3) / dxfij{ri-x(T^j,Vi,ri + {l-x)(Tij,V2;t) 

JQ 

= ^ X! "^'^J' ^-'^ + y ^"^1 /^^^ J da- &{a ■ gi2){a- ■ gi2 fa^ 



/ da; - xcTij, vi,ri + (1 - x) cTij ,V2;t). 



(B19) 



The analysis leading to (|B6p follows from (jB16[) in a similar way with ^pi = miv1/2. However, since energy is not 
conserved in pair collisions the first term on the right side does not vanish. Instead, it represents the collisions energy 
loss w 



{m^v'^ + ■mjv'2 ~ rnivj - m^vl) (ri, Vi,ri + (Tij,V2;t) 

\ X! (1 ^ '^^Aiiicr,^^^ y" c^vi y" dv2 y" d5^e(CT-gi2)(5^-gi2) 



x/y(ri,vi,ri +cr,j,V2;i). (B20) 
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The second term on the right side of (jBlSp gives the coUisional transfer contribution to the flux 



{v'l - vl) / /y(ri - xo-y , vi,ri + (1 - x) (Ty , V2;t) 
Jo 

Vri ■ ^ ^{1 + aij)mi^ijiafj / dvi j dv2 f da e{a ■ gi2){B ■ gu f^ 

■L3 = l J J J 

X [fiji (1 - Qfy ) (a ■ gia) + 2a ■ (^iy vi + MjiVa)] 

X / dx fij{ri- x(T,j,Vi,ri + {1- x)(T,j,V2;t). 
Jo 



(B21) 



This confirms (|B5I) and identifies s'^, which has the equivalent form (obtained by taking half the sum of forms with i 
and j interchanged) 

s'^ = ^ ^ {1 + a,j) m.fijiafj J dvi J dv2 j da <cl{a ■ gi2){B ■ ^12? 
xd- [(1 - ay) {nji - Hij) {a ■ §12) + 4? • (/lyVi + ^J.jiV2)] 

X /y(ri - xcry,vi,ri + (1 - a;)crjj,V2;t). (B22) 

Jo 

APPENDIX C: CHAPMAN-ENSKOG SOLUTION 

As described in the text a normal solution to the kinetic equation is a non-local functional of the hydrodynamic 
fields /i(vi I {yp (t)}). This is equivalent to a function of the fields at a point and all their derivatives at that point 

/,(ri,vi I {yp{t)})^Mvv,{yp{ri,t)};{dr,yp{ri,t) ;..}). (CI) 
If the gradients are small, this function can be expanded in the appropriate dimensionless small parameter 

/.(vi I {yp {t)}) - /r (vi; {yp {ri,t)}) + fP (v^; {yp {r„t)} ; {d^.yp (n, t)}) + • • • (C2) 

where /f'^'' is a function of the fields alone, /l^-* is a function of the fields and linear in their gradients, and so on. 
Thus the kinetic equation can be solved perturbatively by requiring that contributions from common order in this 
gradient expansion vanish. 

To perform this ordering it is necessary to expand the collision operators of (|3.9p 

J„ [ri,vi I /(t)] = afr^ J dv2 J d?e(? • gi2)(? • gi2) 

X - cTy I {n,}) /j(ri,v'/;i)/,(ri - cry,vj,';t) 

-Xij {^1,^1 + o-ij I {"■»}) /i(i"i,vi;t)/j(ri +(T,y,V2;t)]. (C3) 



For the purposes here it is sufficient to go up through first order. The distribution functions evaluated at ri ± cr^ 
become 

/j(vi ; {yp (ri ± cr,^- ,t)};{dr, yp (ri ± 0-^ , i) ;•••})-> (1 ± cr,^- • ) /f"^ (vi ; {yp (ri ,t)}) 
+f^'\^i;{yp {rut)};{^.,yp{r^,t)}) 
= (vi; {yp (ri, 0}) ± (vi; {yp (n, t)})) ■ V,,yp (n, i) 

+/f'(vi;{2/^(ri,i)};R,y/3(ri,i)})- (C4) 

The functional expansion of Xij (ri, ri ± aij \ {ni\) to this order is obtained by a functional expansion of all species 
densities about their values at ri 

Xtj (ri,ri ± cTjj I {n, {t)}) = xlf {(^tj I {n-k (ri,0}) (C5) 
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if 



dr' 



I SXi j (ri,ri ± (Tjj I {rii}) 



s . 
Xif K; {nfc (ri,i)}) + J2 (Vrin^(ri;t)) • J dr' (r' - n) 



^Xij (ri)!"! ±<7-»3 I {n-i}) 
8nfXr',t) 



5n=0 • 



(C6) 



The arrow denotes the leading terms of a Taylor series for (n^(r') — n^(ri)) . The integral can be simplified by noting 
at = the functional integral has translational invariance 



'^X»j(ri,ri±q"ij I {nj}) 
Sne{r',t) 



\sn=o= Fije (ri - r', n ± try - r') 



so 



J dr' (r' 



6xij (ri, ri ± aij \ {ui}) _ /" j„/ /^„/ , 1 \ ^X^j (T^Ty, ±^<T»j | {rii}) 



Sni{r',t) 



\5n: 



-_o= j dr' (^' ± icT, 



5rn{r',t) 



|i5n=0 



(C7) 



(C8) 



1 ainx§^ (ri)}) , [ ^ , , Sxij (T|o-»j, ±|o-»j | {m}) 

2'^'' an,(ri) +y^'"'" 5n,{r',t) 



|(5ra=0 



The expansion for Xij (ri,ri ± cr,j | {m {t)}) becomes 

Xij (ri, ri ± (Tij I {m (i)}) = Xi°^ (o-y-; {nfc (ri, f)}) 

X I (ri) 



1 

1=^ o«^y ■ X^^ri \nne{ri;t) 



i^lnxl"^ (o-»j;{n£ (ri)}) 



(ri) 



+ Iijt{(Jij;{nk (ri,t)}) 



The last line defines Iiji{(Tij; {uk (ri,t)}) as 



Iije{aij;{nk}) = 



'^'^ • " ) 1*"=° 



These results give the expansion of Jy [ri , vi | fi] to first order in the gradients 



(C9) 



(CIO) 



(Cll) 



(C12) 



(C13) 



with the definitions d^X = dX/dr^, 

[Vl I gt, fj] ^ Xi} {<^^fAn^ (ri, t)}) j dv2 1 dae(? • gl2)(? • gl2) 

X [%'5i(V'/)/,(V^') -fli(Vi)/,(V2)] , 



(C14) 



(C15) 
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(ri , V'/ ; (n , V^') + /r (r 1 , Vi ; (r 1 , V2 ) 



(C16) 



Finally, L is a linear operator defined over s dimensional vectors {Xi\ whose components are labelled by the species 



(0) 



(0) 



t(0) 



(C17) 



It remains to choose the magnitude of the external forces to consistently order this expansion. To be specific, 
and for comparison with Ref. [49| it is assumed here to be of first order in the gradients. 



1. Zeroth Order 



At lowest order all gradients of the hydrodynamic fields are neglected, and (|3.8p becomes 

s 

5f /r(vi;fe {v,,t)} = Y,jlf [vi (ri,t)}),/f ({y^(ri,t)}) 



(C18) 



The notation df'^ for the time derivative means that the balance equations are to be used to zeroth order in the 
gradients 



,(0) 



= -C(°)({y/3(ri,t)})r9T/i°^(Fi;{y^(ri,t)}). 
Use of (jCigp in (jClSp gives the zeroth order equation (|6.ip of the text. 



(C19) 



2. First Order 



The kinetic equation for contributions of first order in the gradients is 



3 = 1 



1 * 



vi \nt 



dni 



dy In n£ , 



(C20) 



where Vr = d/dr and Vv = d/dV. The first term on the right side of (jC20[) can be expressed explicitly in terms of 

the gradients, where now d^^^ means that the balance equations are to be used with only terms of first order in the 
gradients 



(9(1) + . Vr,+mriF. . Vv.) /f = mr^F. • Vv + (^../f ) + • V,) y, 

(Vvjf ) • (-p-^Vr,p + Vi • Vr,U) 



(C21) 
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In (IC2T]) . Dl'^' EE dl'-'+U-V and use has been made of the fact that /| ' depends on U only through the combination 
Vi = vi — U, so that 



The pressure gradient can be expressed in terms of the temperature and density gradients 

s 

i=i 

to give 



pT 



ri "-J 



i=i 



Equations (jC20p for the first order distributions, /- , now become 

s 

+ = A,;(Vi;{n,})-Vlnr + ^BnVi;{nJ).Vln, 

If 2 
+Ci^^^ (Vi; {rii}) - d^U,^ + d,^U^ ~ -i^^v^ ■ U 



+A (Vi; {nj) V • U + ^ (Vi; {nj) • F, 



The functions of velocity on the right side of (|C25p are identified as 



(V) = iv.Vv • (v/f )) - ^dyj^^ + \ J: [Vv • (v/f 



1 / , 51nx£ 
- ^ h J , 



f(0) 



(C22) 



(C23) 



(C24) 



(C25) 



(C26) 



(C27) 



(C28) 



1 



A(V) = 3V . Vv/f ^ - ^ ( Ct/ + -^P ) Vv • ( V/„ 



nTd 



f(0) 



(0) 



(C29) 
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Ei(V).-(Vv/rj-(^A.-^J^ (C30) 

Upon deriving (jC26p - ()C30|) . use has been made of the relations 

) ^ -ivv • (v/f )) , du,jr - -dvjf^. (C31) 

The tensor derivative of the flow field d^Un has been expressed in terms of its independent trace and traceless parts, 
using the spherical symmetry of ff'\ e.g. 

+ )V.U, (C32) 

and a similar analysis of the contribution from ^ij,7 ['^ I ^V'/j/j"^]- Equation (|C25p is an inhomogeneous, linear 

integral equation, where the inhomogeneity (the right side) is a linear combination of the the external force and the 
gradients of the hydrodynamic fields. The coefficients of these fields are specified functions of the velocity V. Since 
by definition f^^^ is proportional to the external force and the gradients of the hydrodynamic fields, it must have the 
form 

s 

/f^ ^ A(V)-VlnT + ^Bi(V)-Vlnnj 

+C^nv (V) ^ (^d^Un + d„U^ - ^S^„\/ ■ U 

S 

+I?.(V)V.U + ^5i(V).F,. (C33) 

The unknown functions of the peculiar velocity, A, B:^,Ci^^,,,I?i, and are determined by solving Eq. (jC25p . By 
dimensional analysis, A (V) - v^^fi-'^AUV*), (V) = v^'^i^^'^Bl* (V*), C,,^^ (V) = ^'''^'^^'"''C*^^ (V*), 
A (V) ^ v^'-'^^^^e^'^V* (V*), and Sj (V) = m-^Vg'^'^^^^ fi-'^el* (V*), where £ is an effective mean free path and 
A* (V*), B{* (V*), (V*), V* (V*), and £{* (V*) are dimensionless functions of the reduced velocity V* = V/dq, 
vq — y/2T/m being a thermal speed. Consequently, 



(V) = -C^^^TOtA, (V) = ic^^Vv • (VA (V)) , (C34) 



In addition, 



dl°^Bi (V) = -C^°^T5tB^' (V) = ^C^^^Vv • (VB| (V)) , (C35) 

9f C,,^, (V) = -C^°^TaTC..^, (V) = ic^°^C,,^, + ^C^°Vv • (VC,,^„ (V)) , (C36) 

9f ^A: (V) = -C'^^^^TOtV, (V) = ^C^'^V, + ^C^°^Vv • (VA: (V)) , (C37) 

di'^'^Sl (V) = -(^°^TdT£l (V) = C^"^^- (V) + ^C^^^Vv • (V5^' (V)) . (C38) 



a^°Vlnr = Vaf lnT= -VC^"^ = -ic^"VlnT-^nj^^^Vlnn,. (C39) 
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Since the external force and gradients of the fields are all independent, Eq. (jC25[) can be separated into independent 
equations for the coefficients of each. This leads to the set of linear, inhomogeneous integral equations 



(C40) 



a. 



(C41) 



(C42) 



The linear operator C is 



Notice that (|C40p can be used in (jC4ip to give the equivalent representation for the latter 



(C43) 



(C44) 



(C45) 



(C46) 



This completes the CE solution up through first order in the gradients and first order in the external force. Once 
(16. ip has been solved for ff^ the integral equations for B{, Ci^^n, T^i, and £{ can be solved for fl^\ Then, the 
cooling rate, heat flux, and pressure tensor can be calculated as linear functions of the gradients and the external 
force, and the explicit forms for the transport coefhcients identified. 



APPENDIX D: AN EIGENVALUE PROBLEM FOR £ 

To simplify and interpret the linear integral equations defining the first order solutions it is useful to identify 



a special set of eigenvalues and eigenfunctions for the operator C. Consider the equation for /, 



(0) 



(0) 



ri, Vl I /j , }j 



(Dl) 



Since the temperature occurs through the form (|6.10p . the temperature derivatives can be expressed as velocity 
derivatives 



c(o)vv(v/f)=5:j, 



(0) 
ij 



ri, vi I , 



Noting that C}^^ oc VT, the derivative of this equation with respect to T gives directly 

[CTdrfA =\c,^'''^TdTff'\ 



(D2) 



(D3) 



where use has been made of ([C34| . i.e. Tdrff^'^ = -Vv • (v//°^) /2. An equivalent dimensionless form is 

(rVv • (v/(°)))^ = ic(°)Vv • (v/f ) . (D4) 
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In a similar way, differentiation of (|D2p with respect to each component of the flow velocity gives 

{Cduj^'^)^ = -lc^'^duj'\ (D5) 

Finally, differentiate (|D2p with respect to each of the species densities and noting that the density dependence of 
all quantities occurs only through the /(°)'s and the Xy 's, one gets 

= l{dnX^'' 1/(0, -a„,c^''0^v(v/r) 

= _a„,lnC(°M,(o, ^C('')Vv(v/f). (D6) 

This last form can be simplified by taking into account (|D4p to get 

[C + lnC(°) |^<o, Vv • (v/("))))^ = 0. (D7) 

In summary, there are s + d+1 eigenvalues and eigenfunctions of the operator C. Equations (|D4p and (|D5p identify 
these for the eigenvalue C'"' /2 and the c?-fold degenerate value — C'-'^V^i respectively. Equation (|D7p identifies the 
eigenfunctions for the s- fold degenerate eigenvalue 0. In dimensionless form this eigenvalue problem is written as 

=A('")^'^'"^ (D8) 

The eigenvectors are 

vl/f) = n,9„ jf ) - 2nedn, InC^") 1 <o, *f £ = 1, .., 

= -^Vv • (v/f ) , = -.o5v./r, (D9) 

with the corresponding eigenvalues 

A(™) ^ (o, ..,0, ic(°\ -ic(°\ -^C^"\ -^C^''^) . (DIO) 

These eigenvalues are the same as those of the linearized hydrodynamic equations in the long wavelength limit. This 
provides a direct link between hydrodynamics and the spectrum of the linearized Enskog operator. In addition to this 
physical interpretation, the eigenvalues and eigenfunctions allow a practical formulation of the integral equations, as 
follows. 



AO) AO) 

J i T J i 



1. Biorthogonal set 

Define a scalar product by 

(a,6) = ^ /'dVaI(V)6.(V), (DU) 

1=1 

where the dagger denotes complex conjugation. A biorthogonal basis set is then defined by the eigenfunctions '^'f^^ 
above, and 

V'..= f^,..^,f^^*^-lV^V*V (D12) 
\ni Hi \dm J p J 

where V* = V/uq. The orthonormality condition 

J rfV>„, (V*) ^p, (V*) = Sc.f3 (D13) 



28 



is easily verified. An associated projection operator is given by 



It follows from (|D14p that = V . The corresponding orthogonal projection is 
Consider the quantity (VCX)- 



(D14) 



(D15) 



a j 



(D16) 



Only the projection onto ^'l^^'^^' contributes due to conservation of species number and momentum. It follows then 
that 



(D17) 



The terms with and vanish from symmetry since these all vectors; the last equality follows because C has 

zero trace. Next, note that 



(0) 



(D18) 



or equivalently 



*rc. = *.E / ^-ir^^.[^]-*r^|^E / ^-i^-.-^E^-.t^v^.'/f] 

k 



= '^2^J2 I dviV^2,£, [2?] - P— /C, dv^'f 



(D19) 



This can be used to eliminate the explicit occurrence of the transport coefhcient in the integral equation (jC43p . 
Finally, two additional identities are needed for the proofs of Appendix E 



dv'f, 



(0) p-nT 



P2_LK.,,lv,..(v./f)l, 



(D20) 



(Vv*/r) {vlpY^n.d,,^ (p^nT) = -Vj^ ^^l^ulUdn, + i/.^) (D21) 



APPENDIX E: SOLUBILITY CONDITIONS AND UNIQUENESS 

The results of Appendix D allow proof that the integral equations have solutions and that they are unique. These 
equations have the generic form 



{C-X)X = X, 



(El) 
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where A is one of the eigenvalues (jPlOp . Let relation A" be a solution to (|Eip . Then adding any solution to the 
corresponding homogeneous integral equation also gives a solution 



X' = X + c^, 

where {C — X)'^ = 0. However, the property 

^ - E (V) I dv^„, (V) (V) = 



(E2) 



(E3) 



follows from the fact that the average densities, temperature, and flow velocity are given exactly by the first order 
term fj^^'' , so that contributions to these averages from all higher order terms must vanish. Equivalcntly, (|E3|) implies 



P 



C 
V 



= 0. 



(E4) 



Consequently, the solution to (lEip with the condition (jE4|) is unique. 

To show that solutions exist the integral equations are written in the equivalent form 



(E5) 



dn.i 



(E6) 



(E7) 



(^q{c + ]^C^V^ ^QD,, (E8) 

{Q{c + C)£'),^Qn- (E9) 

These equations are the same as (jC40l) - (jC44[) . The appearance of the factors of Q simply represent a convenient 
rearrangement of those equations, using the identities of Appendix D. They show that the relevant linear operator 
is Q{C — A) where A is one of the eigenvalues (|D10p . The orthogonal projection Q identifies the left eigenfunctions 
with zero eigenvalue as being those of the biorthogonal set il)ai in (jD12[) . According to the Fredholm alternative 
[7l| . solutions to these equations exist if and only if the inhomogeneity is orthogonal to the null space of the left 
eigenfunctions. Here, all the inhomogeneities on the right sides of (|E5[) " (jE9[) appear explicitly orthogonal to this null 
space. Hence, solutions exist and are unique. 



APPENDIX F: DETAILS OF THE CONSTITUTIVE EQUATIONS 



The cooling rate, and fluxes of mass, momentum, and energy are given exactly as explicit integrals of solutions to 
the kinetic equation. Once the CE solution is obtained, approximately to flrst order in the gradients, these expressions 
give the cooling rate and fluxes in the form of the constitutive equations (|5.9p - (|5.12p . The objective of this Appendix 
is to simplify these expressions to the extent possible without making any approximations. This is accomplished in 
most cases by performing solid angle integrations using the results 

j dae{9 •§)(?• g)" = n^'^-'y'-^^g- ^ B^g\ (Fl) 
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/ 



/ 



d?e(?-g)(?-g)"? = B„+ig"i, 



Bn 



da-Q{a • g)(CT • gTdk^i = — —75" {ngkgi + Ski) ■ 

n + a 



J da&{a • g)(CT • g)"CTfcCTfCT„i = g" - 



B. 



n+l 



n + d+l 

where d is the dimension (d > 2), and T (x) is the usual Gamma function 

'1 



YJ>' - 1) dkdedm + dmSki + dkSme + geSkm] , 



r (x + 1) a;r (x) , 



V^, r(i)-i. 



(F2) 
(F3) 
(F4) 

(F5) 



In addition, for the sake of convenience, henceforth we will use the notation gi2 = g. 

To get the collisional transfer contributions to the fluxes, one has to consider the following expansion 

/ dx f ij {ri - x(Tij, vi,ri + {I ~ x) (Tij,-V2;t) = / dxxij {ri - xcTij ,ri + {I - x) (Tij \ {ui}) 
Jo Jo 
x/j(ri - xaij,vi;t)fj{ri + (1 - a;) 0-^, V2;i) 



(0), 



mi 



1 



(F6) 



where Sxij is defined by 



c V^/Y7 r J f^'f' Sxij {ri- x(Tij,ri + {1- x)cr,j \ {n.i}) 

^^'3 = (Vrin£(ri;t)) ■ dx J dr (r - n) 1 5,1=0 (F7) 

The functional derivative is evaluated aA, Sn — and so it depends only on differences of pairs of coordinates, as in 
(|C7|) . A change of variables then makes the dependence on x explicit 



, , _ r ^ ^X,j (n - x(T^J,rl + (1 - x) cr,j | {n,}) 



dr" ( r" + i(Ty - xo-y 



<5n,(r",t) 



i5n=0 



X - - I <T,;i ^ , .. , \Sn=0 + dr r ; , „ \5n=0 



np 



1\ dlnxfMcr.fAnk}) 1 



X - -\nt- 



dui 2 
where Iiji{uij\{nk}) is defined in (jC7p . Finally, then 

1 

-^Xij = ^xi°^ 5Z (^"-1 lnn^(ri; i)) • CTy (/^^^.((Ty ; {n^})) 



(F8) 



(F9) 



1. Cooling rate 



Since C is a scalar, the only gradient contributions are proportional to V-U, and (|4.27p to first order in the gradients 
becomes 



C = C*°^ + Cc/V • u. 



(FIG) 
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with 



i/f (V^)^., • Vv jr(Vi) + 2/f (Vi)I?,(V2) 



(F12) 



Performing the sohd angle integrals gives 



C^"^ - E (1 - 4) / / rfv.<7Vr(V,)/r (V.), (F13) 



2dnT 



2(PnT 



E (1 - 4) / '^vr / rfv. (V.) (g . Vv (VO) 



dnT -^-^ ■' rrii + m,- 

Finally, an integration by parts in the first term of the velocity integrals gives the result quoted in the text 



S f'-^^^^'-''h' /^v,„3,<0.,v,)P,(V,), (F15) 



where the species temperatures are defined by 

= I dvim.FVi°^(V). (F16) 
In the case of mechanically equivalent particles, Eq. (jF15[) reduces to previous results obtained for a monocomponent 

gas til,!?!. 

2. Mass Flux 
The mass fluxes are determined from the definition of (|4.16p 
jo^(ri,i) ^ m, J dvV/f^(ri,v;t) 

= ly dvm,V- (V) Vlnr + ^ (jg^ (V)Vlnnj +5^ (V)F,)j , (F17) 

where the contribution from fj^^"^ vanishes. The transport coefficients according to (jS.lOp are identified as 

A^ = -^ Jd^v-A,{V), (FI8) 
J dvV ■ B{ (V) , (F19) 
j dvV ■ si (V) . (F20) 



mjUjd 



d 
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3. Energy Flux 



The energy flux to first order in the gradients is obtained from Eqs. ()4.18p and (14.19^ as 



q = q 



(F21) 



The kinetic contribution is 



= ^ E/ d^^i^m.V.^Y, ■ A, (Vi) VlnT 

+ ^ E / c^viim^F^Vi . (bI (Vi) V Inn, + (Vi) F, 



(F22) 



The contributions proportional to derivatives of the flow velocity vanish from symmetry. The collisional transfer 
contribution is 

X [- (1 - (^,,- - ^,,) . g) + 4(G,, . a)] [/f Hvi)/f )(V2) + /i')(Vi)/r(V2) 



(F23) 



where Gy = /ly Vi + /ijiV2 . The contribution from in (|F9p vanishes from symmetry. The angular integrals can 
be performed to get 



4^2 



2 + d 



[2(G.,.g).g,+.g2G,,,,] (/f '(Vi)/]^)(V2) + /f ^(Vi)/f (V2)) 



(F24) 



Interchanging the labels i,j and vi,V2 it is seen that the contributions from f^^"^ and /j^-* are the same. For the 
same reason the contributions from dy^^ff^^ and —dy^jf^f^ are the same. The first terms of the integrand give velocity 



the only contributions from dy„fl^^ are those that are scalar functions of the velocities, i.e. those proportional to 



moments of /j^"* of degree one and three, which are proportional to the (partial) mass and kinetic energy fluxes. Finally, 

temperature and species density gradients. The final result is therefore 
" 1 

0== Eid 



, 8B2 



( 2 r^°^ 

2Bi (1 - aij) {n^j - + (d + 2)^ j^^ 

\ nrij ■' mimj •' 



_ (rf + (2m„ - ) + C^V InT + ^ C,:^pV Innp 



(0) 



miirij 



p=i 



(F25) 
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where is defined by Eq. (|F17p and the partial kinetic energy flux is 

q,'= /rfv^F2V/«(V). (F26) 

The constants C,^ and Cf^^ are 

Cl^^^Jd^.J rfv2 [- (1 - a.,) ifi,, - ^i,,)g' + 45 (g • G,,)] /f ^ (VOTOt/]"' (V2), (F27) 

CL = ^ j / ^^^2 [- (1 - a„ ) (m., - H^) 5' + 4g (g • Gy )] /i°HVi)np9„^/f (V2). (F28) 
The expression of Cfj can be simplified when one takes into account the relation 

Tdrff^ ( V) = - i V V • ( V/j°) (V)) , (F29) 
and integrates by parts in (jF27p . The result is 

Cl = j dvi j dv2/f (Vi)/f (V2){5G?^.+5-^(g-G,,f + (1 + M,.).9(g-G,,) 

+H^t^tj9^ + ^ (1 - ) (a^j* - ) [9 is ■ ) + 3^] I ■ (F30) 

On the other hand, no significant further simplification of Eq. (|F28p is possible until f^^'' is specified in detail. 
The heat fiux is seen to have the form (|5.11[) . 

s 

q(r,<) -> -AVT- (^'-D^.^.^lnn, +L,,F,) , (F31) 

so the transport coefficients now can be identified 

A = A*" + A'', Dq^ij = D^'jj + D'^j^j, Lij = + L.^-. (F32) 

The kinetic parts are, 

f '^v^y^V . A (V) , (F33) 

i=l i=l 

= "d^ / ^v^V^^V . Bl (V) , (F34) 

L%^-\j d.'^V^M-el{Y), (F35) 
while the coUisional transport parts are given by Eqs. (|7.14p - (|7.16p . 

4. Momentum Flux 

The momentum fiux to first order in the gradients is obtained from (j4.20p - (|4.22p 

P^x^p!^^ + P'^^, (F36) 
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where 



(F37) 



X [/f (Vi)/f (V2) + 2/f)(Vi)/ji)(V2) 



(F38) 



A contribution to (|F38p proportional to the density gradients from the expansion of Xij(ri — xcr^j, ri + (1 — x) cTij) 
vanishes from symmetry. For similar reasons, the only gradients contributing to both (|F37[) and f IF38[) are those from 
the flow field. The terms proportional to Vi in (IC33P also do not contribute due to the orthogonality condition (jE4|. 
The solid angle integrations can be performed with the results 



-> S^xnT +lflj d^im.,VixVi^C^,0,, (Vi) (^dpU,, + d^Up - • , 



(F39) 



(0) 71(0) ^ 



ijXij nifij 



mi rrii 



B2 

'd + 2 



■ My (1 + «y) x!?«.4 / dv2m,V2^V2xC,.p^ (V2) (dpU^ + d^Up ~ ^Sp^^W ■ U 
i,j=i ^ 



1 S, 



23 + d 



An integration by parts in the velocity integral, and use of fluid symmetry gives, finally 



(0) rr,{0) - 



mi m," 



d + "y) "iCT^'^i / dv2mjV2yV2\Cjj:Sf, (V2) (dpUf, + df,Up - ^^/J^V • U 

+ E (1 + "^. ) 4^4 / / rfv2/f )(VO/r (V2)5 



2d (d + 2) ^ 



aAC/7 + d^Ux - ^S^xV ■ v) + ^-J^'J^aV • U 



(F40) 



(F41) 



The pressure tensor therefore has the form (|5.12p . and the pressure, shear viscosity, and bulk viscosity are identified 
in terms of their kinetic and coUisional transfer contributions. Their expressions are given by Eqs. (|7.18p ~ (l7.22p . 
respectively. 
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